//=============================================================================
//
//   Copyright (C) 2020 All rights reserved.
//
//   Filename:		mod.rs
//
//   Author: Wang Zhecheng - wangzhecheng@yeah.net
//
//   Date: 2020-09-24 01:03
//
//   Last Editors: Wang Zhecheng - wangzhecheng@yeah.net
//
//   Last modified:	2020-11-06 13:20
//
//   Description:
//
//============================================================================
pub mod duncanchang;
pub mod elastic;

use crate::arith::tensor::{Strain, Stress};
use crate::dimension::DimType;
use crate::precision::{Real, TOL};
use crate::util::{record::Record, HasName, Indexable};
use approx::relative_eq;
use nalgebra::{allocator::Allocator, DefaultAllocator, MatrixN};
use serde::{Deserialize, Serialize};
use std::marker::PhantomData;
use std::rc::{Rc, Weak};

pub fn build_material<DT: DimType>(data: &serde_json::Value) -> Rc<dyn MaterialLike<DT>> {
    let name = data["type"]
        .as_str()
        .expect("Failed to get Material's type name")
        .to_string()
        .to_lowercase();
    log::trace!("New material {} has been built.", name);
    match name.as_str() {
        "elastic" => {
            let mat: elastic::Elastic<DT> = serde_json::from_value(data.clone())
                .expect("Failed to deserialize Elastic Material");
            Rc::new(mat) //as Rc<dyn MaterialLike<DT>>
        }
        "duncanchang" => {
            let mat: duncanchang::DuncanChang<DT> = serde_json::from_value(data.clone())
                .expect("Failed to deserialize DuncanChang Material");
            Rc::new(mat) //as Rc<dyn MaterialLike<DT>>
        }
        _ => panic!("{} is not a valid material type", name),
    }
}

pub trait MaterialLike<DT: DimType>:
    HasName + Indexable + MaterialCommonParameter + std::fmt::Debug
{
    ///   get non-linear dee matrxi according to given element state parameter refers to [elastic dee matrix] in elastic material or [non-linear elastic dee matrix] in hypo-elastic material or [elasto-plastic dee matrix] in elasto-plastic material1
    fn get_dee_matrix(&self, s: &Stress<DT>, rec: &Record<DT>) -> MatrixN<Real, DT::DimVoigt>
    where
        DefaultAllocator: Allocator<Real, DT::DimVoigt, DT::DimVoigt>;
    /// update stress history related parameters for given sigma and epsilon (only override in hardening/softening material for changing yield function)
    fn update_paremeter(&self, s: &Stress<DT>, dpe: &Strain<DT>, rec: &Record<DT>);
    /// get the value of the yield function
    fn get_yield_function(&self, s: &Stress<DT>, rec: &Record<DT>) -> Real;
    /// get [D_{elastic}] matrix according to Linear E and v
    #[inline]
    fn get_elastic_matrix(&self) -> MatrixN<Real, DT::DimVoigt>
    where
        DefaultAllocator: Allocator<Real, DT::DimVoigt, DT::DimVoigt>,
    {
        let e = self.get_elastic_modulus();
        let v = self.get_elastic_poisson();
        DT::get_elastic_matrix(e, v)
    }
    /// integrate sigma through epsilon path with sub-stepping method
    ///
    /// Reference:
    /// @article{ doi:10.1002/nme.1620240505, author = {Sloan, S. W.},
    /// title = {Substepping schemes for the numerical integration of elastoplastic stress–strain relations},
    /// journal = {International Journal for Numerical Methods in Engineering},
    /// volume = {24},
    /// number = {5},
    /// pages = {893-911},
    /// doi = {10.1002/nme.1620240505},
    /// url = {https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.1620240505},
    /// eprint = {https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.1620240505},
    /// abstract = {Abstract This paper describes two substepping schemes for integrating elastoplastic stress–strain relations. The schemes are designed for use in finite element plasticity calculations and solve for the stress increments assuming that the strain increments are known.  Both methods are applicable to a general type of constitutive law and control the error in the integration process by adjusting the size of each substep automatically. The first method is based on the well-known modified Euler scheme, whereas the second technique employs a high order Runge–Kutta formula. The procedures outlined do not require any form of stress correction to prevent drift from the yield surface. Their utility is illustrated by analysis of typical boundary value problems.},
    /// year = {1987}
    /// }
    #[inline]
    fn integral_through_epsilon_path(
        &self,
        de: &Strain<DT>,
        s: &mut Stress<DT>,
        rec: &mut Record<DT>,
    ) where
        DefaultAllocator: Allocator<Real, DT::Dim, DT::Dim>
            + Allocator<Real, DT::DimVoigt>
            + Allocator<Real, DT::DimVoigt, DT::DimVoigt>,
    {
        // 1. trial delta sigma
        let ds = Stress::<DT>::from_voigt(self.get_elastic_matrix() * de.to_voigt());
        // 2. if trial sigma lies in elastic zone
        if self.get_yield_function(&(*s + ds), rec) <= 0 as Real {
            *s += ds;
        }
        // 3. if trial sigma lies in plastic zone
        else {
            // 3.1 get alpha: the ratio of elastic/all in trial delta sigma
            let alpha = self.get_ratio_of_nonlinear(&s, &ds, &rec);
            // dpe: delta plastic epsilon
            let dpe = (1 as Real - alpha) * *de;
            // 3.2 add delta elastic sigma. find the elastic sigma at yield surface.
            *s += ds * alpha;
            // 3.3 iteration to find final sigma/epsilon
            let mut d = 0 as Real;
            let mut dt = 1 as Real;
            while d < 1 as Real {
                dt = if dt < (1 as Real - d) {
                    dt
                } else {
                    1 as Real - d
                };
                // ddpe: trial detal plastic epsilon
                let ddpe = dt * dpe;
                d += dt;
                // 3.3.1. first trial
                let ds1 = Stress::<DT>::from_voigt(self.get_dee_matrix(&s, &rec) * ddpe.to_voigt());
                // 3.3.2. get new sigma and state paramenter (refers to new yield function)
                let s1 = *s + ds1;
                let rec1 = rec.clone();
                self.update_paremeter(&s1, &ddpe, &rec1);
                // 3.3.3. second trial
                let ds2 =
                    Stress::<DT>::from_voigt(self.get_dee_matrix(&s1, &rec1) * ddpe.to_voigt());
                // 3.3.4. average
                let ds_mid = (ds1 + ds2) * 0.5;
                // 3.3.5. convergence check
                let diff = (ds2 - ds1) * 0.5;
                let relative_diff = diff.value().norm() / (*s + ds_mid).value().norm();
                if relative_diff < TOL {
                    d += dt;
                    // add plastic delta sigma
                    *s += ds_mid;
                    self.update_paremeter(&s, &ddpe, &rec);
                } else {
                    let mut beta = (TOL / relative_diff).sqrt() * 0.8 as Real;
                    beta = if beta > 0.1 as Real {
                        beta
                    } else {
                        0.1 as Real
                    };
                    beta = if beta < 2.0 as Real { beta } else { 2 as Real };
                    dt *= beta;
                }
            }
        }
    }
    /// get ratio of non-linear part to a given increment stress
    /// Reference:
    /// @ misc{ports2001finite,
    /// title={Finite element analysis in geotechnical engineering},
    /// author={Ports, DM and Zdravkovic, L},
    /// year={2001},
    /// publisher={Thomas Telford, London}
    /// }
    #[inline]
    fn get_ratio_of_nonlinear(&self, s: &Stress<DT>, ds: &Stress<DT>, rec: &Record<DT>) -> Real
    where
        DefaultAllocator: Allocator<Real, DT::Dim, DT::Dim>,
    {
        let mut alpha_last: Real = 0 as Real;
        let mut alpha_curr: Real = 1 as Real;
        let mut alpha_next: Real;

        // secant iteration
        while relative_eq!(alpha_last, alpha_curr, epsilon = TOL, max_relative = TOL) {
            let f_curr = self.get_yield_function(&(*s + alpha_curr * *ds), rec);
            let f_last = self.get_yield_function(&(*s + alpha_last * *ds), rec);
            if relative_eq!(f_last, f_curr, epsilon = TOL, max_relative = TOL) {
                alpha_next = 0 as Real;
            } else {
                alpha_next = alpha_curr - f_curr / (f_curr - f_last) * (alpha_curr - alpha_last);
            }
            alpha_last = alpha_curr;
            alpha_curr = alpha_next;
        }
        alpha_curr
    }
}

#[derive(Debug, serde::Deserialize)]
#[serde(untagged)]
pub enum Material<DT: DimType> {
    Elastic(elastic::Elastic<DT>),
    DuncanChang(duncanchang::DuncanChang<DT>),
}

macro_rules! material_impl {
    (CommonParameters:
        $($cparam:ident$(($crename:expr))?:$cparamtype:ty),*;

    $($mat:ident$(:
        $($sparam:ident$(($srename:expr))?:$sparamtype:ty),*)?);*) => {

        #[derive(Debug, Clone, Serialize, Deserialize)]
        $(pub struct $mat<DT:DimType>{
            #[serde(skip_serializing_if = "Option::is_none")]
            id: Option<usize>,
            #[serde(skip)]
            _marker: PhantomData<DT>,
        })*
    };
}

material_impl!(
    CommonParameters:
        volume_weight: Real,
        friction_angle: Real,
        dilation_angle: Real,
        cohesion: Real,
        perm("permeability"): nalgebra::Vector3<Real>,
        zeta: Real,
        k_ic("K_IC"): Real,
        elastic_modulus: Real,
        elastic_poisson: Real;

    Elastic;
    DuncanChang:
        rf("Rf"):Real,
        k("K"):Real,
        n:Real,
        g("G"):Real,
        f("F"):Real,
        d("D"):Real
);

pub trait MaterialCommonParameter {
    fn get_volume_weight(&self) -> Real;
    fn get_friction_angle(&self) -> Real;
    fn get_dilation_angle(&self) -> Real;
    fn get_cohesion(&self) -> Real;
    fn get_permeability(&self) -> nalgebra::Vector3<Real>;
    fn get_zeta(&self) -> Real;
    fn get_k_ic(&self) -> Real;
    fn get_elastic_modulus(&self) -> Real;
    fn get_elastic_poisson(&self) -> Real;
}

// pub trait MaterialPoolLike<DT: DimType>: std::fmt::Debug {
//     fn append_material(&mut self,rhs)
// }
pub trait Linked2MaterialPool<DT: DimType> {}

///
/// has_material trait
///
pub trait HasMaterial<DT: DimType> {
    fn get_material(&self) -> Rc<dyn MaterialLike<DT>>;
    fn set_material(&mut self, rhs: Weak<dyn MaterialLike<DT>>);
}

#[test]
fn serde() {
    use crate::dimension::PlaneStrain;
    let data = r#"
        {
            "id": 2,
            "type": "elastic",
            "volume_weight": 19,
            "friction_angle": 30,
            "dilation_angle": 30,
            "cohesion": 10,
            "permeability": [1e-5,1e-5, 1e-5 ],
            "zeta": 0.1,
            "K_IC": 2.9,
            "elastic_modulus": 1000000.0,
            "elastic_poisson": 0.3
        }"#;
    let mat: Material<PlaneStrain> = serde_json::from_str(data).unwrap();
    println!("{:?}", mat);
    // assert_eq!(mat.get_volume_weight(), 19 as Real);
    // assert_eq!(mat.get_elastic_poisson(), 0.3 as Real);
}
